
rm(list=ls(all=TRUE))

library(Zelig)
library(miscTools)

set.seed(123)

setwd("replication_archive/Cellphones/")


### Load the data
load("data_public.rda")
data <- new
rm(new)


### Run original models from Table 1: Logit (m2) and Fixed Effects OLS (m6)
data <- subset(data, select=c("conf2008_dum", "cell_dum_07", "pre2000_count", "bdist1", "capdist", "pop2005", "mnt", "irri", "gcppc00", "cell_dum_07", "cow"))
data$cow <- as.factor(data$cow)
data <- data[complete.cases(data)==T,]


### Run the sensitivity analysis
etavec <- seq(0, 0.022, length.out=50)

results2 <- results6 <- NULL

for(i in 1:length(etavec)){
	res2 <- res6 <- NULL
	for(j in 1:500){
		# create data given eta
		data$conf2008_true <- NA
		data$conf2008_true[data$conf2008_dum==0 & data$cell_dum_07==0] <- rbinom(length(data$conf2008_true[data$conf2008_dum==0 & data$cell_dum_07==0]), 1, etavec[i])
		data$conf2008_true[data$conf2008_dum!=0 | data$cell_dum_07!=0] <- data$conf2008_dum[data$conf2008_dum!=0 | data$cell_dum_07!=0]

		# run the models with the simulated data
		m2 <- zelig(conf2008_true ~ cell_dum_07 + pre2000_count + bdist1 + capdist + pop2005 + mnt + irri + gcppc00, model="logit", data=data, cite=F)
		m6 <- zelig(conf2008_true ~ cell_dum_07 + pre2000_count + bdist1 + capdist + pop2005 + mnt + irri + gcppc00 + factor(cow), model="ls", data=data, cite=F)

		# save the coefficients and standard errors, as well as the number of units with cellphone coverage that have a violent event
		est2 <- m2$getcoef()[[1]][2]	
		err2 <- sqrt(m2$getvcov()[[1]][2,2])
		res2 <- rbind(res2, c(est2, err2, sum(data$conf2008_true)))

		est6 <- m6$getcoef()[[1]][2]	
		err6 <- sqrt(m6$getvcov()[[1]][2,2])
		res6 <- rbind(res6, c(est6, err6, sum(data$conf2008_true)))				

	}
	results2 <- rbind(results2, c(etavec[i], colMedians(res2)))
	results6 <- rbind(results6, c(etavec[i], colMedians(res6)))
}


colnames(results2) <- colnames(results6) <- c("phimax", "pointest", "stderr", "number")


### create confidence intervals
results2 <- as.data.frame(results2)
results6 <- as.data.frame(results6)

results2$lowerci95 <- results2$pointest + qnorm(0.025)*results2$stderr
results2$upperci95 <- results2$pointest + qnorm(0.975)*results2$stderr

results6$lowerci95 <- results6$pointest + qnorm(0.025)*results6$stderr
results6$upperci95 <- results6$pointest + qnorm(0.975)*results6$stderr



### FIGURE 20

### Panel (a)
results <- results2

quartz(type="pdf", width=5, height=5, file="output/cellphones_m1.pdf")
par(mar = c(4,4,0.3,0.3), mgp=c(2.5,1,0), family="CMU Serif")
plot(1:length(results$pointest), results$pointest, type="n", ylim=c(min(results$lowerci95), max(results$upperci95)), xlab = "Percent Additional Events in Non-Cellphone Units", ylab="Coefficient", xaxt="n")
polygon(c(1:length(results$pointest), rev(1:length(results$pointest))), c(results$lowerci95,rev(results$upperci95)), col="grey", border=NA)

points(1:length(results$pointest), results$pointest, type="l", lwd=3)
axis(1, at=seq(1, 50, length.out=5), labels = seq(0, 80, length.out=5), las=2)
abline(h=0, col = "black", lwd=2)
lines(c(1,1), c(results$lowerci95[1], results$upperci95[1]), lwd=3)
points(1, results$pointest[1], pch=16)
dev.off()


### Panel (b)
results <- results6

quartz(type="pdf", width=5, height=5, file="output/cellphones_m5.pdf")
par(mar = c(4,4,0.3,0.3), mgp=c(2.5,1,0), family="CMU Serif")
plot(1:length(results$pointest), results$pointest, type="n", ylim=c(min(results$lowerci95), max(results$upperci95)), xlab = "Percent Additional Events in Non-Cellphone Units", ylab="Coefficient", xaxt="n")
polygon(c(1:length(results$pointest), rev(1:length(results$pointest))), c(results$lowerci95,rev(results$upperci95)), col="grey80", border=NA)

points(1:length(results$pointest), results$pointest, type="l", lwd=3)
axis(1, at=seq(1, 50, length.out=5), labels = seq(0, 80, length.out=5), las=2)
abline(h=0, col = "black", lwd=2)
lines(c(1,1), c(results$lowerci95[1], results$upperci95[1]), lwd=3)
points(1, results$pointest[1], pch=16)
dev.off()



### FIGURE 19
results <- results2

# total events where C=0
sumnull <- sum(data$conf2008_dum[data$cell_dum_07==0])

# total events where C=1
sumone <- sum(data$conf2008_dum[data$cell_dum_07==1])

# total number of cells where C=0
cellnull <- length(data$conf2008_dum[data$cell_dum_07==0])


quartz(type="pdf", width=5, height=5, file="output/cellphones_descr.pdf")
par(mar = c(4,4,0.3,0.3), mgp=c(2.5,1,0), family="CMU Serif")
plot(etavec, (results$number-sumone)/cellnull, type="l" , ylim=c(0, (0.002+max((results$number-sumone)/cellnull))), lwd=3, xlab=expression(eta), ylab=c("P(T'=1)"))
dev.off()



